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MULTIGRID METHODS FOR DIFFERENTIAL 
EQUATIONS WITH HIGHLY OSCILLATORY 

COEFFICIENTS* 

Bjorn Engquist Erding Luo 
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SUMMARY 

New coarse grid multigrid operators for problems with highly oscillatory coefficients are 
developed. These types of operators are necessary when the characters of the differential 
equations on coarser grids or longer wavelengths are different from that on the fine grid. 
Elliptic problems for composite materials and different classes of hyperbolic problems are 
practical examples. 

The new coarse grid operators can be constructed directly based on the homogenized 
differential operators or hierarchally computed from the finest grid. Convergence analysis 
based on the homogenization theory is given for elliptic problems with periodic coefficients 
and some hyperbolic problems. These are classes of equations for which there exists a 
fairly complete theory for the interaction between shorter and longer wavelengths in the 
problems. Numerical examples are presented. 

INTRODUCTION 

Multigrid methods are usually not so effective when applied to problems for which the 
standard coarse grid operators have significantly different properties from those of the fine 
grid operators [1,3,7-9,11-12]. In some of these problems the coarse grid operators should 
be constructed based on other principles than just simple restriction from the finest grid. 
Elliptic and parabolic equations with strongly variable coefficients and some hyperbolic 
equations are such problems. One feature in these problems is that the smallest eigenvalues 
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do not correspond to very smooth eigenfunctions. It is thus not easy to represent these 
eigenfunctions of the coarser grids. 

We shall investigate elliptic equations with highly oscillatory coefficients, 

^ 0 j-( x )^ u ‘(*) = -ft®)* = aj ( x ’ 7^ 

with cij(x,y) strictly positive, continuous and 1-periodic in y. This is one class of the 
problems discussed above for which there exists a fairly complete analytic theory such that 
a rigorous treatment is possible. This homogenization theory describes the dependence of 
the large scale features in the solutions on the smaller scales in the coefficients [2,11]. We 
shall consider model problems but there are also important practical applications of these 
equations in the study of elasticity and heat conduction for composite materials. 

In this paper we analyse the convergence of multigrid methods for equation (1) by 
introducing new coarse grid operators, based on local or global homogenized forms of the 
equation. We consider only two level multigrid methods. For full multigrid or with more 
general coefficients the homogenized operator can be numerically calculated from the finer 
grids based on local solution of the so called cell problem [2]. 

In a number of numerical tests we compare the convergence rate for different choices of 
parameter and coarse grid operators applied to a two dimensional elliptic model problem. 

The convergence rate is also analyzed theoretically for a one dimensional problem. 
If, for example, the oscillatory coefficient is replaced by its average, the direct estimate 
for multigrid convergence rate is not asymptotically better than just using the damped 
Jacobi smoothing operator. The homogenized coefficient reduces the number of smoothing 
operations from 0(h~ 2 ) to 0(/i“ 10/7 log h). When h/t belongs to the set of Diophantine 
numbers [4], ergodic mixing improves the estimate to 0(h~ 6 / 5 log h ). The step size is h and 
the wave length in the oscillating coefficient e. 

These results carry over to some but not all hyperbolic problems. A numerical study 
of using hyperbolic time stepping with multigrid in order to compute a steady state gives 
similar results to the elliptic case. 

TWO DIMENSIONAL ELLIPTIC PROBLEMS 

Elliptic problems on the form (1) will be considered, 

- V • a'(j:,y)Vu £ = f(x,y), (x, y) G fi = [0, 1] x [0, 1], (2) 

subject to Dirichlet boundary condition u e | an = 0. The function a e (x,y ) = a(x/e,y/e) is 


176 


strictly positive and 1-periodic in x and y. From homogenization theory [2] follows, 


max I u, — u| — * 0, as e — » 0. 

where u satisfies the following effective equation, 

subject to the same boundary condition. Here, 

r d /c ■ 

— J , <S 2 ) 1 2 } ^ > J 

and the periodic functions tz 3 are given by, 

da^Si , s 2 ) 


( 3 ) 


-V 5 -a(s 1 ,s 2 )V s Ac J = 




j = 1,2. 


For the numerical examples we shall choose a special case with diagonal oscillatory coeffi- 
cient, 

* 7 * 7 / 

( 4 ) 


1-1/ 

a f (x,y) = a( - ). 


From (3), we know that the corresponding homogenized equation is, 

(y + a)d 2 u d 2 u (fi + a) d 2 u 

+ n — izrr. = J{ x ’y)i 


2 d 2 x 


dxdy 


2 d 2 y 


( 5 ) 


where fi = m(l/a f ) _1 and a = m(a c ). Here, the mean value m(f) of a e — periodic function 
is defined as, 

m(f) = - f f{x)dx. 

6 JO 

For convenience, we introduce a brief notation of a N x N block tridiagonal matrix T , 


T — Tridiag[T i _ 1 ,T l ,T i+1 }^ NxN ^ 


T u 

T 2 \ l' 22 t 23 


Tnn-i Tj 


NN 


111 


Numerical Algorithm 

The discretization of (2) combined with (4) is 


- D* a h D x u h . - Dy_b h D y u h . = f h .. 

+ *3 - tj + - V J % 3 


(6) 


where a*. = -\~ Vj), = a £ (x,- - y 3 + \ ), ij = 0, • • • , N. D + and D _ are forward 

and backward divided differences, respectively; h = denotes the grid size. Using vector 
notation, we can rewrite (6) as 


L (lh U i>h 




where 


( 7 ) 


L Clh = JjT ridiag[B^_ v A .J, 

A) = + a£ + 6* + b h {j _ v -ay (7V _ 1)x(W _ 1) 

Bj is a diagonial matrix, denoted by B*> — Diag[—b^ .]^ N _ x) X (iV— i) an ^ 

= ( U n> U 21’ ' ' " > U jV-ll’ ' ’ ’ ’ U 1N-V U 2N-1' 1 ' ’ U N-1N-1 ) 

F _( fh fh ... fh . . . f, h fh . . . fh )T 

* t,h V J 1 1 ’ J 21 ’ ’■'JV-11’ ’ J\N-V J 2N-V ’-'TV-lW-D 

For simplicity, we only consider the two-grid method. Denote the full iteration operator of 

this method by M. It is defined by, 


M = S-’(I-It l LJlFL„ h )S\ 


( 8 ) 


where the restriction and interpolation operators are given, as denoted below, by the weight- 
ing restriction and bilinear interpolation operators, respectively, 


I 1* = — 
h 16 
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2 
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r h 


r 1 

r — * | ^ 
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The smoothing iteration operator S is based on the damped Jacobi iteration, 

S = I — u>h 2 L tll . 

The coarse grid operators Ljj is one of the following operators: 


( 9 ) 
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Global Homogenized operator: which is the discretized form of (5) 


—0.5 (/x + a)D* + £Fuij + (/x - d)D*D y u l3 - 0.5(// + a)D» = / tJ . 
Written in matrix form, 

^ P > -Sf ] ( f - i ) X ( f _ 1 ) > 


( 10 ) 


where , _ . _ 

/i + a + a 

= Tridiagl — ,2{g + a), — J ( ^_ 1)x( ^_ 1) , 

rr , r a — a a + a a — a, 

Bf = Tndiag{ — , —£—]<$ -i)x(f -ip 

Local Homogenized operator: L^ has the same form as (10), except the entries for BR 
coming from the local discretized values of a t {x — t/), 


Lit = -j^Tridiag[Bf_^M\ S"] ( ^_ 1)x( ^_ 1} , 


(11) 


where 


with 


A f = Tridiagl-aH^a?^. + a” + + MJ, -«£](*_, )x( f _ 1} , 

B f = Tridiag\—c¥_ v — 6^, c^ r ] ( ^_ 1)x( ^_ 1) , 


m 

u 


«? 17-1 ' *J-1' 


a^. + a^. t + 26(ab a h . a — a 

13 7 ~u -i' // _ tj ij-l V h __ i 

’ a «i “ 4 ’ * " 2 ’ 


%j,c 2 ) is defined to be . 

Reduced Load Homogerihed operator^ L H has the same form as in (11), except here we ig- 
nore the cross term D*Dq. That means BH is a diagonal matrix, BH = Diag[-bV ] ( £_i) x (£_! 


Lji = ~Tridiag[BR_ v A^Bf] ( E_ 1)x( K_ x) (12) 

Sampling operator: L c H has the exact form as L c h , but values a-, b tJ are defined on the 
coarse grids, 

L h = L CiH (13) 



Variational operator: 


(14) 


J H 




Numerical Results 


In practice, it is not always easy to calculate the spectral radius p(M). Therefore, we 
study the mean rate [14] of convergence under different coarse grid operators L H . The 
mean rate of convergence is defined by 


P = ( 


\\ L C,h Ui ~ AjU xJly 

WL^hU 1 - AH* 


(15) 


where i is the smallest integer satisfying \\L t h u' — fh.\\h < 1 x 10 6 . 

In Figure 1, a c (x — y) — 2. \ + 2 sin(27r(i — y)/e). We plot p defined by (15) as a function 
of 7 by taking e = y/2h, and u in (9) is 0.095. 






Figure 1 : p as a function of 7 . Dotted line is for (10), solid line for (11), dashed line for (13), and dashdot for 
( 12 ), + for (14). ( 1 . 1 )-( 1 . 4 ) are for different number of grid points N . 
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It is clear that the coarse grid operators derived from the homogenized forms (10) 
and (11) are superior. The effect is more pronounced for large 7 when the eigenspace 
corresponding to large eigenvalues of L t h is essentially eliminated. For the practical low 
7 case, a study of the impact of the choices of I fa and Iff is needed. In this paper we are 
concentrating on the asymptotic behavior (large 7 ). Different I ^ and Iff operators are 
briefly discussed for the one-dimensional problem. 

In Figure 2, we plot p as a function of the variable a, where L H (a ) comes from the 
discretized operator —a^D x +D x _ + a(p — a)D x D q — b^-D\D v _^ <x> = 0.095 and e = y/2h. 



Figure 2: p as a function of a. Here TV = 64 and 7 = 12. denote p under the different choice of Normal 
and Local Homogenized coarse grid operator, respectively. 

From Figure 2, we get further evidence of the importance of using the correct homoge- 
nized operator. Techniques based on one-dimensional analysis does not contain the mixed 
derivative term [ 1 ]. 

In order to isolate the influence of the coarse grid approximation we have kept the 
smoothing operator fixed. It obviously also affects the performance. If we use Gauss Seidel 
iteration method in (9), the convergence rate can be improved. In Table 1 , we test the 
same coefficient a e (x — y) = 2.1 + 2 sin(27r(x — y)/ e). Taking N = 128, e = y/2 h, we compare 
the convergence rate by choosing damped Jacobi iteration and Gauss Seidel iteration. 


7 

5 

6 

7 

8 

9 

10 

11 

12 

Jacobi 

.5929 

.5519 

.5173 

.4863 

.4579 

.4349 

.4140 

.3950 

G-S 

.4545 

.4221 

.3922 

.3703 

.3491 

.3304 

.3158 

.3008 


Table 1 : Spectral radius, two dimensional case 
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ONE DIMENSIONAL PROBLEMS 


The one dimensional equation is useful as a model for which a more complete mathe- 
matical analysis is possible. 

Consider the following one-dimensional elliptic boundary value problem with a periodic 
oscillatory coefficient, 

d du, 

— a c (x)— — = 1, 0 < x < 1, u £ (0) = u e (l) = 0, (16) 

dx dx 

where a e (x) = a(^) and satisfies the same assumption as above. As e — > 0, u t converges 
strongly in L to the solution u of the homogenized equation, 

d 2 u 

— a— j = 0 < x < 1, a = m(l/a c ) L (17) 

dx 1 

Subject to the boundary conditions <^(0) = </>(l) — 0. 

Numerical Algorithm 

Let the difference approximation of (16) be of the form: 

- a £ (x J+ i)(u^ +1 - uj) + a ‘(*,-_p(uJ - uj_,) = 1, 1 (18) 

In matrix form, (18) can be written as 

L ( ,hU h = 1, = {u\ ■ ■ w^_ 1 ) T 

where r 

L t ,h — T^Tridiagl—a^, a t +x](Ar-i)x(N-i) (19) 

with a.j = a c (xj — |). 

We first consider a two-grid method by applying standard restriction, standard inter- 
polation operators and Jacobi smoothing iteration. 

The coarse grid operator L H will be one of the following: 

Homogenized operator: 

m(l/a ‘) -1 . , 

= Tndiag[— 1 , 2 , — l](^-i) X (^-i) (20) 

Averaged operator: 

l h = ^rTridiag[- 1,2, -l] ( £_ 1)x( £_ 1) (21) 
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Sampling operator: has the exact form as L c but only every second a 3 value is used, 


L h — L tH (22) 

Variational operator: 

" (23) 


Convergence Theory 

The theorem below on the convergence rate is too pessimistic in the number of smooth- 
ing iterations necessary. However, the analysis still gives insight into the convergence 
process and the role of homogenization. With Lr replaced by averaging (21) the same 
analysis results in 7 = 0(h~ 2 ) which means that multigird does not improve the rate of 
convergence over just using Jacobi iterations. This follows from the effect of the oscillations 
on the lower eigenmodes. It should also be noticed that in the case (ii), the solution of L t h 
is much closer to those of Lr , see [ 11 ]. 

Theorem 1 Let M be defined as in (8), with L H defined by (11). There exists a constant 
C such that , 

p(M) < p 0 < 1 , 

when either one of the following conditions is satisfied : 

(%) 7 > Ch-'-VUnh 

(ii) the ratio of h to e belongs to the set of Diaphantine number } and 7 > Ch~ l ~ l / S lnh. 

For details of the proof, see [ 10 ]. An outline is as follows . Separate the complete eigenspace 
of L t h into two orthogonal subspaces, the space of low eigenmodes and that of high eigen- 
modes. After several Jacobi smoothing iterations in the fine grid level, the high eigenmodes 
of the error are reduced, and only the low eigenmodes are left. Combining eigenvalue analy- 
sis with homogenization theory [ 11 ], one may realize that the low eigenmodes of the original 
discrete operator are close to those of the corresponding homogenized operator. We then 
approximate them by the corresponding homogenized eigenmodes and correct these in the 
coarse grid level. 
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Numerical Results 


In Figure 3.1 and Figure 3.2, a £ (x) = 2.1 + 2 sin( 27 r x/e). We plot the analogous graph 
to Figure 1. Here e = V2h and u in (9) is 0.1829. In Figure 3.3 and Figure 3.4, a c (x) = 
2.1 + 2 sin( 27 ra;/e + tt/ 4). Here e = 4/i and u? = 0.1585. 
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Figure 3 : p as a function of 7 . Dotted line is for (20), solid line for (21), dashed line for ( 22 ), and dashdot for 
(23). ( 3 . 1 )-( 3 . 4 ) are for different number grid points N . 

In Figure 4, with the assumptions in Figure 3. 3-3. 4, we plot a c (x) and the approximation 
of (18) under the choices of coefficients in Figure 3. 



Figure 4: (4.1) and (4.2)are the graphs for a‘(x), where * are the discretized values. (4.2) is the solution. 
Dashed line is for (17). Dashdot line is for “m(a < )ti rx = 1 and line with circles is for (18). 
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In Figure 5, we plot p as a function of the variable a, where L H = aA ff . In (5.1), 
a e (x) = 2.1 + 2sin(2'K x j e) , u = 0.1829 and e = \f2h\ In (5.2), a c (x) = 20.1 ifx/t — [x/e] G 
(0.7, 0.9); otherwise, 0.1, to = 0.0373 and e = 4/i. 




Figure 5: p as a function of a. The homogenized value ah = m(l/a e ) 1 and the arithmetic value av = m(a e ) 
are given. Here 7 = 10 and N = 256. 


In Figure 6, we present the convergence u t — » u, as e — ► 0 by giving the numerical solutions 
of (16) and (17). Recall that our goal is to solve the oscillatory problem and to use the 
homogenized operator only for the coarse grids. 




Figure 6: Solid lines are the approximations for (18), dashed lines are the solutions for (17), respectively, 
when c = 0.2 in (6.1) and e = 0.1 in (6.2). Here N = 500. 





HYPERBOLIC PROBLEMS 


Time evolution of a hyperbolic differential equation can be used for steady state computa- 
tions. This is common in computational fluid dynamics, [6]. In multigrid this means that 
hyperbolic timestepping replaces the smoothing step. There are fundamental differences 
with standard multigrid for elliptic problems but some of our earlier discussions carry over 
to the hyperbolic case. The dissipative mechanisms for hyperbolic problems are mainly the 
boundary conditions. Consider using the model problem, 


d 2 u t 

1 W 


d , ,du t 


0 < x < 1 


(24) 


cis the smoothing equation in multigrid for the numerical solution of (16), subject to the 
boundary conditions 

u e (0) = 0, = 0. (25) 


dx 


The equation (24) must have boundary conditions which are dissipative but reduce to (25) 
at steady state, see [5], 


m«(o, 0 = o, 


du t (l,t) 

dt 


+ \! a "{ x ) 


duJXt) 

dx 


= 0 . 


(26) 


The initial condition should support the transport of the residual to the dissipative 
boundary x = 1, 

u e (x,t) = u°(x) given, 


t (x, t + Af) = u°(x) — A t^Ja<^x)D*u°(x). 


Note that the initial condition approximates the transport equation u t + ^/au x = 0. The 
difference approximation of (24) needs a low level of numerical dissipation. 

The homogenization theory of [2] is also valid for equation of the type (24). A numerical 
indication is seen in Figure 7. 

The positive effect of multigrid on the convergence rate does not carry ever to problems 
for which the steady state is hyperbolic or contains hyperbolic components. If, 


du t 

dt 


du, „du. 

+ a— 1 - = 0 

dx dy 


is used for the equation, 


u c is 


du, „du, . . 

a—+P~r- = $, x, y £ [0, 1], 

1 — periodic in y, u t (0,y,t) — a e (y). 
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The coarse grid operator must resolve all scales of a c to required accuracy in order to 
produce multigrid speed up. More on this phenomena will be reported elsewhere. 

Numerical Results 

In Figure 7, take 50 smoothing steps. Coefficient a(x/e) is the same as in Figure 1. 






Figure 7: (7.1) Solutions: Solid line is the solution of steady state; Dashed line for homogenized solution; 
Dashed dot line for average solution. (7.2) Residue as function of two level multigrid cycles. (7.3) Approximate 
solutions after each two level cycle. (7.4) Approximate solutions for time evolution equation. 


CONCLUSION 

Elliptic equations and some hyperbolic equations with highly oscillatory coefficients 
have been studied. We have shown that the homogenized form of the equations are very 
useful in the design of coarse grid operators for multigrid. 
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The evidence is from a sequence of numerical examples with strongly variable coefficients 
and to some extent from theoretical analysis. The result is clear in the asymptotic regime 
of many smoothing iterations. 

The impact on the coarse grid operator from the numerical truncation error and the 
interpolation operator needs to be asessed in order to improve the performance in the 
regime of very few smoothing iterations per cycle. 
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